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Abstract 

We have developped a procedure for the search of periodic signals in the 
data of gravitational wave detectors. We report here the analysis of one year of 
data from the resonant detector Explorer, searching for pulsars located in the 
Galactic Center (GC). No signals with amplitude greater than h = 2.9 10 -24 , 
in the range 921.32-921.38 Hz, were observed using data collected over a time 
period of 95.7 days, for a source located at a = 17.70 ± 0.01 hours and 5 = 
— 29.00±0.05 degrees. Our procedure can be extended for any assumed position 
in the sky and for a more general all-sky search, even with a frequency correction 
at the source due to the spin-down and Doppler effects. 



1 



1 Introduction 



Periodic or almost periodic gravitational waves (g.w.) are emitted by various 
astrophysical sources. They carry important information on their sources (e.g., 
spinning neutron stars, accreting neutron stars in binary systems) and also on 
fundamental physics, since their nature can test the model of General Relativity 
||. ||. The main feature of continuous signals which allows them to be detected 
is that, despite the weakness of the signal (compared to typical amplitudes for 
bursts), it is possible to implement procedures that build up the signal to noise 
ratio (SNR) in time. The natural strategy for searching for monochromatic 
waves is to look for the most significant peaks in the spectrum. In this case 
the SNR increases with the observation time t f, s . In fact, as t Q b s increases, 
the frequency resolution of the spectrum also increases - the frequency bin gets 

smaller bv = thus the content of the noise power in each bin decreases 

tobs 

with t bs, while the signal is not dependent on the length of observation timeQ 
For a periodic signal of amplitude h at the frequency v the squared modulus 
of the Fourier Transform provides h 2 with a noise contribution of 2 Sh(v) 8v, 
being Ship) the two-sided noise power spectrum of the detector (measured in 
Hz- 1 ). 

Thus the SNR for periodic signals is: 

Eq.|l| holds if the instantaneous frequency of the continuous signal at the 
detector is known. The analysis procedure in this case is "coherent", since 
the phase information contained in the data is used and the sensitivity (in 
amplitude) increases with the square-root of the time. However in some cases 
it may be impossible, for various reasons (see later in sect. 3.1), to perform a 
single Fourier Transform over all the data. This means that the observation 
time has to be divided in M subperiods, such that the spectral resolution of the 
spectra becomes bv = M/t Q b s and the corresponding SNR is M times smaller 
than that given by eq. [j]. 

The M spectra can be combined together by incoherent summation, that is 
by averaging the square modulus. In this case the final spectral resolution is 
again bv but there is still some gain as this procedure reduces the variance of 
the observed noise in the bin. We obtain: 

SNR' = - tob ° (2) 
2S h (v)VM 



1 Th.e noise power density, i.e. the power spectrum, is independent of t b s , while the signal energy 
density grows with t b s 
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In general, if the signal is monochromatic but frequency modulated due 
to the detector-source relative motion, processing techniques exist which can 
recover the sinusoidal case if the source direction is known. One of the standard 
ways of detecting such signals is through appropriate resampling of data, better 
known in the radio astronomy comunity (where this technique is commonly 
used) as "data stretching" (see for example ||). In the case of radio pulsar 
searches, the location of the source is usually known (the data come from a 
radio telescope pointing in a particular direction) but some parameters of the 
system need to be estimated and this is done by a "timing solution which is 
phase coherent over the whole data set"||. 

However for gravitational waves, especially when searching a large param- 
eter space, it is doubtful that strategies developed for radio pulsar searches 
can simply be adapted: the expected low SNR values for g.w. signals really 
modifies the nature of the search strategies that can be employed. In recent 
ms pulsar searches, for example in (3|, the signal is strong enough to allow 
suspected pulsars to be identified by visual inspection of the results of the final 
stages of the analysis procedure. 

For gravitational waves the study of the implementation of optimum analysis 
procedures is still in progress @, |5|, 0, §1 • 

In the present analysis, we analyzed the data of the Explorer detector to 
search for ms pulsars located in the GC, assuming their intrinsic frequency to 
be constant over the analysis time. 

The procedure we used in this study relies on a data base of FFTs, computed 
from short stretches of data (short with reference to the effects of the Doppler 
shift, as will be described later in this paper). These short FFTs are then 
properly combined together to provide a new set of FFTs with higher frequency 
resolution, representing the signal in the frequency range selected for the study. 

The paper is organized as follows: in section ^ we briefly review the charac- 
teristics of the detector during the 1991 run, in section || we describe the main 
aspects of the procedure; in section |] we present the results obtained. The Ap- 
pendices clarify some aspects of the analysis procedure and discuss the extent 
to which the constraints that we have introduced to describe the procedure can 
be relaxed in order to account for different sources. 



2 The Explorer detector 

The Explorer detector is a cryogenic resonant g.w. antenna located at CERN, 
at longitude 6° 12' E and latitude 46°27' N. The apparatus and the experimental 
set up of the antenna during the 1991 run have been described in || and some 



results of the data analysis for burst detection are given in [1C]. 

The system has two resonance frequencies {v_ = 904.7 Hz and u + = 
921.3 Hz in 1991) where the sensitivity is highest. Fig.[l| shows the behaviour 
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of the two resonance frequencies with time. 

Fig.^ shows the hourly averages of the energy sensitivity (SNR=1) to mil- 
lisecond bursts, expressed as effective temperature T e ff in kelvin, obtained with 
an adaptive Wiener filter 0. The relation between T e j f and the amplitude of a 
ms burst is [[| h = 8 10~ 18 \/T e fj (T e ff in kelvin). 

For periodic waves the sensitivity of a bar detector at the above resonances 



is given by Jl^, 13]: 



= 2M . 1Q - 25 T 2300 ^900^1^, 

V 0.05 K M Q i/ t obs v ; 

where T is the bar temperature, M its mass, Q the merit factor, vq the res- 
onance frequency of the mode and t b s the observation time. After one year 
of effective observation, the minimum detectable h (amplitude detectable with 
SNR=1), using the nominal parameters of the Explorer detector (T = 2 K, 
M = 2300 kg, Q = 10 6 ), is 

h = 2 • 10~ 25 (4) 
in a bandwidth of ~ 2 Hz around the two resonance frequencies and 

h ~ 2 • 10" 24 

in a bandwidth of 16 Hz between the two resonances. For the NAUTILUS 
[14j or AURIGA ju| detectors (with T = 0.1 K,Q = 10 7 ) we get a value 



ft, ~ 1.5 • 10 26 at the resonances. 



3 Main features of the analysis procedure 

In the search for continuous signals there are a number of issues that need to 
be kept in mind regarding the signals that might be present, the apparatus and 
the quality of the data. As far as the source is concerned, it is not possible to 
set up a single procedure capable of searching over all types of periodic signals. 
In this analysis, we concentrated on periodic signals such as those expected 
from isolated neutron stars with weak spin down, i.e. we ignored the spin down 
parameters ^. Moreover we did not consider the effects of proper accelerations 
of the source. Thus our model assumes that the frequency behaviour of the 
signal exclusively depends on the Doppler effect caused by the motion of our 
Earth-based detector relative to the source location. Let us quote some basic 
figures: 



2 the sensitivity obtained with a matched filter was, on average, better by a factor of two. The 
comparison between the two filtering procedures is shown in [|1J . 

3 however it is possible, using our procedure, to take into account the spin down. This will be the 
next step in the development of the procedure. 
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Figure 1: Top: lower resonance (minus mode) frequency against time (in days of 
the year). Bottom: upper resonance (plus mode) frequency against time. A slow loss 
in the electrostatic charge of the transducer is probably responsible for the frequency 
drift 
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Figure 2: Hourly averages of the Explorer sensitivity to millisecond bursts, expressed 
as noise temperature (K) as a function of time (in days of the year) 
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the Doppler effect has two periodic components (see Appendix A). The first 
one, due to the revolution motion of the Earth over a period of 1 year, produces 
a maximum spread of 

Au rev = 0.1986 V cos Ped mHz (5) 

and a maximum time derivative 0.197 • 10 -10 V cos (3 ec i Hzs -1 where V is the 
intrinsic frequency of the signal and (3 ec i is the ecliptical latitude of the source. 
The second one, due to the rotation of the Earth over a period of 1 sidereal 
day, produces a maximum spread of 

A^Voi = 0.00308 V cos Pter cos 5 mHz (6) 

and a maximum time derivative 1.12 • 10~ 10 V cos ftter cos 5 Hzs -1 . fiter is the 
terrestrial latitude of the detector and 5 the declination of the source. 

The observation is also affected by modulation in the amplitude. This is due 
to the varying orientation of the detector with respect to the source because 
of the Earth's motion. It may also be a consequence of the polarization of the 
wave. 

As shown, for example, in |l6|], this modulation spreads the signal power 
across side bands, spaced at 1/24 hours. The amplitude modulation observed 
using a resonant bar detector is given by the geometrical part of the detector 



cross section 17 



E = S • <S>(6, e, <t>) = S • sirfO ■ (— ^— + ecos 2 (2(j))) (7) 

where e is the percentage of the wave linear polarization (e=l when there is 
linear polarization, e=0 when there is no polarization), 

IT C C 

is the (two-sided) mechanical part of the cross section, (M=bar mass, v= sound 
velocity in the bar), 9 is the angle between the bar axis and the wave direction 
of propagation, 4> is the angle between the bar axis and the wave polarization 
plane. The cross section is maximum when the parameters are 9 = ir/2,(fi = 
,e = 1. If the source location and the polarization state are known, it is 
possible to demodulate the amplitude of the observed signal. 

A major consideration in developing the analysis is that the operation of 
the detector is not continuous and the noise is not stationary. An example of 
this is given in Fig.|| which shows power spectra computed over two hours of 
data in November and September 1991. There are several lines from periodic 
disturbances which are not stationary, and the noise level differs between the 
two spectra. 
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Figure 3: The figure shows two power spectra of the detector, both obtained during 
periods of "good" operation of the apparatus: during November (left) and September 
(right) 1991. The y-axis is Sh ■ 10 40 in units of 1/Hz. The x-axis is the frequency in 
[Hz]. Comparing the two power spectra, it is easy to see the non stationarity of the 
system 

A final remark needs to be made about the time accuracy required to per- 
form a coherent search over a given time period. A rough calculation, shown 
in Appendix B, shows that, in a coherent search, the global timing error must 
be less than 10 -4 — 10~ 5 s/year. 

3.1 The length of the FFTs in the data base 

Our frequency domain data base consists of "elementary spectra", each ob- 
tained by performing the FFT (Fast Fourier Transform) of a given number of 
samples of the data, over a duration to, recorded by our detector. The signal 
sensitivity of each spectrum, according to eq. |], depends on the duration to- As 
our observations are affected by the Doppler shift, we have chosen the duration 
to of the "elementary spectra" to be the longest possible compatible with the 
requirement that the signal should "appear as monochromatic" during to- This 
is expressed as: 



*o opt = max to such that 




Clearly some assumptions about the expected z> must be made. In principle, 
in order to achieve a higher SNR in the short spectrum, some preprocessing 
could take place by setting a coarse grid on the parameter space (i.e., the part 
of the sky that is being investigated) and performing suitable data stretching 
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for each point in that parameter space. In this way a signal coming from that 
parameter space would appear as monochromatic in the resulting spectrum. As 
a consequence the dimension of the data base is increased by a factor equal to 
the number of points in parameter space, but there is a gain in SNR because 
of the higher spectral resolution. 

As stated above, we restricted our analysis to the case where the only fre- 
quency changes are due to the Doppler effect of the detector motion relative to 
the source. 

As shown in Appendix A, the time duration must be 



with 

a = V ■ 11.244 • 1CT 11 cos (3 ter Hz (daily rotation), 

b = V ■ 1.977 • 10~ n cos (5 ec i Hz s _1 (annual orbital motion), 

and V is the source intrinsic frequency. 

Thus, to construct the elementary spectra of our data base, we chose a 
duration of to = 2382.4 s=39.7 minutes, corresponding to 2N = 131072 sam- 
ples, recorded with sampling time of 18.176 ms. With this choice of to, the 
maximum Doppler frequency variation, Su^ = 0.28 mHz, is smaller than the 
resulting frequency bin bv = 0.419 mHz. 



The header of each elementary spectrum of the database contains various infor- 
mation about the original data. This allows stretches of data that are noisier 
than others to be vetoed or weighted differently and thus best exploits the 
potential of the data. 

Some of the information contained in the header relates to the data struc- 
ture, some of it to the operational status of the detector and some to data 
quality. For example, the date and time of the first sample of the data series 
that the FFT is computed from is stored, along with the frequency resolution 
of the FFT and the type of time-domain windowing used. There are also sys- 
tem parameters that vary in time: the frequencies of the two modes and of 
the calibration signal, the level of brownian noise and merit factors at the two 
modes, the wide band noise level and the status of the operation flags (normal 
operation, working around the antenna, liquid helium refilling). Some of this 
information was used to set a threshold for vetoing data (at 40 minutes, which 
is very reasonable due to typical non stationarities of the detector) 





3.2 The FFT header 
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3.3 The procedure for combining spectra 

For the targeted search described here the basic FFTs are combined coherently 
to improve the final sensitivity. The following is an outline of how this is done 
(details are given in Appendix C). 

• take an FFT. Let bv be the frequency resolution, 2N the number of sam- 
ples and B the bandwidth of the detector; 

• take the data from n' bins in the frequency band Au of the actual search; 
n' = NAu/B; 

• build a complex vector that has the following structure: 

- the first datum equal to zero 

- the next n' data equal to those from the selected bins of the FFT 

- zeroes from bins n' + 1 up to the nearest subsequent bin numbered with 
a power of two. Let us say that this way we have n bins. 

- zeroes in the next n bins 

So, we end up with a vector that is 2n long. 

• take the inverse FFT of the vector. This is a complex time series that 
is the "analytical signal" representation of the signal in the band Au. It 
is shifted towards lower frequencies and it is sampled at a sampling rate 
lower by a factor 2N/2n compared to the original time data 0. The time 
of the first sample here is exactly the same as the first datum used for the 
data base and the total duration is also that of the original time stretch. 
There are fewer data because the sampling time is longer. 

• repeat the steps outlined above for all the M FFTs; 

• if they all come from contiguous time stretches simply append one after 
the other in chronological order. If they are not all contiguous set to zero 
those stretches where data are missing. 

• To correct for the Doppler effect Q from a given source, multiply each 
sample of the sequence by 

eX pi-mi)) (io) 

U are the times of the samples and <j>{ti) = Au £>(t)dt. 

Auj£>(t) is the Doppler correction, in angular frequency, at the time t 

of the i th sample: Aui£,{t) = — ui s , where cj£>(i) is the frequency 



4 the construction of the analytical signal is a standard procedure of lowpass filtering for a 
bandpass process. In fact the analytic signal is zero on the left frequency plane, thus avoiding 
aliasing effects in the lowpass sampling operation [Q 

5 we could also take into account other causes of frequency shifts such as those affecting the 
intrinsic frequency of the source. 
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observed at the detector, due to the Doppler effect from a given source 
that emits at a constant frequency u s . t s is the start time of the FFT 
being constructed. 

We note that the frequency correction is performed on the subsampled 
data set, and this is one of the advantages of the procedure. 

• Perform the FFT of the 2n ■ M data thus obtained. 

• Finally, take the squared modulus of the FFT thus obtained. This is 
the power spectrum of the original time series, in the band Av, with the 
full spectral resolution 5v/M. A signal exhibiting frequency variability 
smaller than the variability we have corrected for, should appear wholly 
within a single frequency bin, and its resulting SNR will be that of Eq.|| 

3.3.1 An example of the procedure for combining spectra 

The procedure was tested on simulated signals added to the data. We shall now 
briefly review the results of these tests. Such simulations, although simple in 
principle, present practical design problems which demand extreme care in the 
implementation. The simulated signal is constructed in the time domain and 
then it is handled in exactly the same way as the real detector data (details on 
the data handling procedures are given in the Appendices). Each FFT of the 
signal is then added to the corresponding FFT in the data base. 

s(nAt) = h(nAt) sin (<f)(nAt) + O ) (H) 

where 

At is the sampling time, n = 0, 1... .131072, <pQ =initial phase, 

rnAt 

4>{nAt) = / oj D (t)dt (12) 



where LUr>(t) is the frequency at the detector due to the Doppler shift at time 



t. Using the discrete form of eq.12 we may write the phase at time tf. 

(pi = (f>i-i + u Di At 

We report here an example in the absence of noise (we set the detector 
FFTs to zero, before adding them to the simulated signal). 

FigJ3] shows the comparison of the two power spectra obtained from a source 
assumed to be in the GC, emitting at 921.3 Hz, before and after Doppler 
removal. It is clear that the spread and the shift in the signal frequency (top 
figure) have been properly corrected (bottom figure). Here the observation time 
is 36 hours and the frequency resolution is 6A...fiHz. 

The level of the signal, after Doppler removal, is that which would be ex- 
pected (h = 1.0 / \J Hz). An accurate analysis of the residual error after Doppler 
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Figure 4: Simulation of a signal (at 921.3 Hz) from the GC, over 36 hours of data. 
Top: Shiy) of the simulated data. Bottom: S^iy) after Doppler removal. The y-axis 
is Sh • 10 40 in units of 1/Hz. The x-axis is the frequency in [Hz]. 

removal (this error is defined as the instantaneous difference between the time 
signal after correction and the time signal in the absence of modulation) showed 
that this residual error was always less than 0.7%. 
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4 The analysis of the Explorer data 



We analyzed data taken during the period between March and December 1991. 
These data sets comprise 4954 FFTs from March to July and 4384 from August 
to December. After a preliminary analysis of the features of the spectra, with 
particular reference to their sensitivity performance, we decided to veto the 
spectra with brownian noise larger than 7.8 K (i.e. three times greater than 
the expected value of 2.6 K ||), With this criterium we vetoed 807 spectra, 
that is ~ 10 % of the total. 

A comment on the accuracy of the timing of the data is necessary at this 
point: the absolute time recording had an indetermination of the order of 
10 — 20 ms ||] at the beginning of each new run. This was due to the fact 
that, although the time was checked against the Swiss time signal HBG with an 
accuracy of a few ms, the software procedure at the start of each run introduced 
an imprecision of ~ 10 — 20 ms. 

On the other hand we are confident about the precision of the Rubidium 
clock, which was used to determine the sampling time. As a consequence we 
could combine coherently only data obtained from one single acquisition run. 

The strategy for the analysis procedure was thus the following: 

• choose the frequency bandwidths to be analyzed, and calculate new - 
higher resolution- FFTs on each new run, for each of these chosen band- 
widths 

• choose the coordinates of the source location and correct for Doppler effect 
and for amplitude modulation, using the procedure described in section 
|3[ We ignored possible polarization of the waves. 

This analysis was focused on possible sources in the GC, at a = 17.7 hours, 
S = —29.0 degrees. 

To calculate the Doppler shift we used the JPL ephemerides (JPLEPH.405) 
and software routines from the US Naval Observatory (NOVAS). The amplitude 
modulation was removed from the data by matched filtering, multiplying the 
data by the factor sin 4 6{t) (in other words, the data were weighted on the basis 
of the source-detector direction) . 

In 1991 we collected data over 51 separate runs and therefore, applying 
the procedure outlined above, we obtained 51 separate FFTs. Each one has a 
different frequency resolution, according to its length in time (due to the timing 
error between different runs). The analysis of the 51 FFTs could only be done 
by combining their information "incoherently", thus combining the spectra. 
Obviously, this reduced the sensitivity of the final analysis (see Eq. |2|). 

We report here the results of the analysis of the data around the frequency 
of the plus mode (the mode in Fig.[l], bottom). 
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4.1 The analysis of 95.7 days of data 



First of all we give an example using data over one week in June 1991. The 
observation time is t b s = 7.05 days from day 159.8 (June, 8th). The analyzed 
data are in the bandwidth 921.32-921.38 Hz. The level of the noise is (1.2 ± 
0.7) 10 -24 , in good agreement with that expected for Explorer (using eq.|l], with 
T = 2.6 K and Q = 10 6 , we get 1.9 • 10~ 24 ) 

The Doppler correction needed for signals from the GC was applied to the 
data. The highest peak is h = 5.2 10 -24 . 

In the case that no Doppler correction was applied to the data there is a 
very high peak, at the level h = 1.2 10~ 23 , which disappears when the Doppler 
correction based on a source in the GC is applied. It is most likely that the peak 
was due to the apparatus. The GC Doppler correction spreads this contribution 
over several frequency bins. 

We started the analysis using only the data from three consecutive runs in 
May and June. 

We averaged the corresponding spectra over a total observation time of 
t bs = 21.177 days from day 128.53 (May, 8th). The analysis was done in the 
frequency range 921.32-921.38 Hz, where the antenna noise spectral density was 
flat, as shown in Fig.|5[ The figure is normalized in terms of the amplitude h 
that would give SNR = 1 for sources in the GC. The h level is (1.6±0.5) 10~ . 
We notice that the average level for the 21.177-day period is roughly (within 
the observed non stationarities) the same as that for the 7.05-day period, and 
the standard deviation decreases, as it should do, by a factor of the order of 



= \/3 (the small difference is due to the non stationarities during the 



three time periods). 

No spectral lines were detected with amplitude (at the detector) greater 
than h = 4.1 10" 24 . 

To set an upper limit on the amplitude of possible signals from GC in the 
chosen bandwidth, we decided to check the efficiency of detection, given the 
noise of the detector. We therefore added signals with different amplitudes and 
phases to the data, using data without Doppler correction, since the efficiency 
of detection, on the average, does not depend on the Doppler effect. 

We added four different families of signals, each family consisting of 20 
sinusoids with the same amplitude but different phases. 

Fig.^ shows the histograms. The histograms report the data and the four 
different families of simulated signals. On the x-axis we have values ranging 
from h = 1.0 10~ 24 to 1.6 10~ 23 , with an interval of 0.1 10~ 24 . 

The nominal amplitudes of the added signals are 1.4 • 10 -23 , 8.7 • 10~ 24 , 
5.8 • 10 -24 and 3.0 • 10~ 24 . If we consider only the three families of higher 
signals, the histograms show very clearly that all these signals have been well 
detected. Even the smallest of these 60 signals is well above the standard 
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deviation of the data (the smallest three signals have h = 4.4 10 -24 , which is 
roughly 10 times the noise standard deviation). Thus, for these signals, on a 
time basis of only 21 days, the efficiency of detection is 1. On the contrary, the 
histograms show that the efficiency of detection for the signals at the lowest 
level (3.0 • 10~ 24 ) is very poor. 

Thus, on the basis of 21 days of data, we exclude the possibility that, in the 
GC, there are sources having a spindown age r > 3 • 10 7 years emitting signals 
with frequency in the range 921.32-921.38 Hz and strength (on Earth) greater 
than or equal to h = 5.8 10~ . 

For the eleven longest runs between May and December the observation 
times range from 7.7 to 12.8 days, giving a total effective observation time of 
95.72 days. We can average these spectra, after adding the necessary zeroes 
to obtain the same virtual resolution (this produces a change in the SNR and 
therefore a re-calibration of the spectra is needed). 

The start times of these eleven runs, are days 128.53, 137.29, 159.82, 171.43, 
213.94, 225.32, 301.61, 312.37, 323.46, 339.68. Fig.0 shows the average h for 
the frequency range 921.32-921.38. In the frequency bandwidth 921.32-921.38 
Hz over these 95.7 days the noise level is (1.2 ± 0.2) • 10~ 24 , well in agreement 
with the expected value. No lines with amplitude greater than h = 2.9 • 10 -24 
are apparent. 

The standard deviation is a factor 2.5 lower than the standard deviation 
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Histogram of the data and of the added signals {Tobs=21 days) 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

h X10" 3 

Figure 6: Histograms in the bandwidth 921.32-921.38 Hz. t f, s = 21 days. The 
histograms report the data (magenta) and the four different families of simulated 
signals (violet ,yellow,green,red). The y-axis numbers above five have not been plotted. 
The x-axis ranges from 1.0 10 -24 to 1.6 10~ 23 , with interval 0.1 10 -24 . The histograms 
show clearly that the efficiency of detection for the simulated signals corresponding 
to the yellow, green and red plots is 1. 
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Observation time=95.7 days 
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Figure 7: Average h from GC obtained using data from May to December (95.7 days) 
in the bandwidth 921.32-921.38 Hz. The x-axis is the frequency, in Hz-921 Hz 

obtained using only 21 days, thus we expect the efficiency of detection over the 
95 days be of the order of the unity even for signals of h = 2 — 3 10~ . 

Thus, we exclude the possibility that, in the GC, there are sources hav- 
ing a spindown age r > 10 8 years, emitting signals with frequency in the 
range 921.32-921.38 Hz and strength (on Earth) greater than or equal to h = 
2.9 10" 24 . 

4.2 Point registration on the spectra 

The analyzed period consists of 51 runs, leading to 51 spectra of different 
resolution. It is not convenient to average these spectra as done before for the 
eleven longest ones, because now their duration are very different one from each 
other. They can be analyzed using various methods- for example, by looking 
for patterns in the time evolution of their spectral lines. However this kind of 
analysis would require algorithms which are rather more involved than those 
used in the present analysis (such algorithms are presently under investigation 
[||, The analysis here is also complicated by the fact that the different 
spectra have different resolutions and thus different SNRs, for any given signal. 

We restricted our search to a source in the GC emitting a signal at constant 
frequency during the observation time. We have tracked all the local maxima 
in each spectrum obtained by setting a threshold^ at SNR ~ 4. If a spectral 
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line from the GC were present, it should show up in all the spectra (at various 
SNRs) at the same frequency. 

Fig. H shows (top) the time- frequency plot of the selected maxima and their 
histogram (bottom). The resulting histogram is flat and hence no evidence of 
straight horizontal lines is present in the data. However, as shown in FigJ^, the 
sensitivity of this analysis is much poorer than the previous method, as almost 
all these selected maxima have amplitude greater than 1CP 22 . 
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Figure 8: The incoherent analysis over 51 spectra from May to December. The 
upper plot shows the time-frequency behaviour of the local maxima selected in each 
spectrum. The lower plot is their histogram versus frequency. 
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Histogram of the detected amplitude of the maxima 




log h 10**(-20) 

Figure 9: The figure shows the histogram of the (log) amplitudes of the selected local 
maxima in the 51 spectra 
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5 Conclusions 



A first analysis of the data obtained with the Explorer detector in 1991 was 
performed, with the aim of searching for continuous g.w.. The analysis was 
limited to the frequency range 921.32-921.38 Hz, which contained the plus res- 
onance of the detector, where sensitivity was highest. Doppler corrections on 
the g.w. frequency were made under the assumption that the source was still 
in the GC without any intrinsic frequency spin down. 

No signals were observed with amplitude greater than h = 2.9 10 -24 , using 
data collected over 95.7 days, for a source located at a = 17.70 ± 0.01 hours 
and 5 = —29.00 ± 0.05 degrees, having a spindown parameter r > 10 8 years 
(that isp < 1.7- 10~ 19 s/s). 

The procedures adopted here can be applied to any assumed position in the 
sky of a g.w. source, for a greater frequency range, or even for a frequency 
correction at the source due to spin down and Doppler effects. 
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APPENDIX A: Frequency resolution of the 
basic FFTs 



Here we want to show how we chose the spectral resolution for our spectral 
data base [16, 17|. 

If V is the emission frequency of the source then the observed frequency is 
given by the (simplified) formula 

v = F(l ± (y/c)) (13) 

where v = \v r ot + is the relative velocity of the antenna and source, due to 
v ro t, daily rotation of the Earth, and to v orb , orbital motion around the Sun. 
Orbital motion, annual: 

Vorb = \Vorb\ ' COs(Q orb t) (14) 

\v or b\ = tt orb R = 2irR /T orb = 2vri? /(365 • 86400) = 29.772 km/s (15) 

R = 149504201 km; n orb = 2vr/(86400 • 365) = 1.99 • l^rad/s; \v orb \/c ~ 
0.9933 • 10" 4 . 

_ = „(l/c) cos^^- (16) 

= Ro^lrb - cos(tl orb t) = v orb Vt orb (17) 



^JT-lmax = V (\v orb \/c) COS/3 ed orb = V (\v orb \/c) COS Peel — (18) 

dt l orb S 

Hz 

~ V ■ 1.9778 • 10~ n cos (5 ed (19) 

s 

Daily rotation: 

Vrot = \v T ot\ ' COs(Q rot t) (20) 

broil = ttrotRr = 2irR T /T rot = 2irR T /864O0 = 0.4636 km/s (21) 
R T = 6378.388 km; Q rot = 2vr/86400 = 7.2-1Q- 5 rad/s; \v rot \/c ~ 1.5467-1Q- 6 . 



^JT^lmax = V (\v rot \/c) COS (3fer ^rot = V (\v rot \/c) COS Pter — ( 22 ) 

at l TOt s 

Hz 

~V- 11.244 • 10" 11 COS Pter 

s 
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We will use here the symbols \Au max \ rot and e \Av max \ or b for the maximum 
frequency shift, in the time interval AT. 



I Au max \ ro t 

(P AT) 
if AT less than 1 day. 



11.244 • 10" 11 cos I3 ter s- 1 (23) 



AiW l°^ = 1.977 -10- 11 cos p ecl s- 1 (24) 



(V AT) 
if AT less than 1 year. 

(\Au max /AT\ rot + \Au max /AT\ orb ) ■ T (T < 1 day) (25) 
(|Azw/ATU + |Aiw/ATU) • T < (26) 
l/T >T-{a + b) 

where 

a = V ■ 11.244 • 1CT 11 cos/3 ter (daily rotation) 

b = V ■ 1.977 • 1CT 11 cos /3 ec ; (annual orbital motion) 

T 2 < l/(o + 6) 
T < 8.6972 • 10 4 /VF s 
If we consider F = 1000 i^z 

T < 2749.9 s (27) 

In our data, the sampling time is St f = 18.176 ms then, using T = 2382.35 s, 
that is 0.6617 hours, we have a frequency resolution 

S Vf = 0.41975 mHz 

that is 2 17 = 131072 samples in each periodogram. 



APPENDIX B: Timing precision 

We note here that precision in both the absolute time and the sampling fre- 
quency during the observation time i {, s is crucial when combining the basic 
FFTs to construct one very long spectrum. 
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If v is the frequency of the wave, t = \ the wave period, e v = the relative 
frequency error (Av represents here the wave frequency indetermination due to 
the acquisition and analysis procedures) then 

At = e u ■ t 

After time t m we have a phase error of 

— = t m Av = t m ve v 

Z7T 

Hence, if we want a phase error ^ <~ 0.01 in t m = 1 year = 3.1 • 10 7 s we 
need e u <~ 3- 10~ 10 - 
Hence, if v = 1 kHz: 

e u <~ 3 • 10~ 13 
At <~ 10~ 5 s/?/ear 

APPENDIX C: Practical issues in the con- 
struction of the data base 

• Each FFT is computed using 2N data, sampled with sampling time At. 
The data are windowed, in the time domain, before the Fourier trans- 
form. This means that the data yi are multiplied by Wi = A — B cos (i) + 
C cos (2i), where i = (0, 2N — 1) • 2tt/(2N — 1). In the present analysis we 
have used a Hamming window, that is A = 0.54, B = 0.46 and (7 = 0. 

• The FFTs are stored in units of strain/ \/Hz, and are normalized so that 
their squared modulus is the spectrum. 

• The basic FFTs of the data base overlap for half their length. The time 
duration of each FFT is to = 2NAt, and a new FFT is done after time 
to/2. This is important since it avoids distortions in the final time domain 
sequence - this is the well known "overlap-add" method, described in many 
data analysis textbooks. For the Explorer detector we thus have 110 FFTs 
over 36 hours with 5v = 0.41... mHz. 

• We select the frequency bandwidth to be analyzed and we add zeroes to 
construct the analytical signal. These data should be to a power of two, 
to allow use of a Fast Fourier algorithm. The chosen bandwidth should 
be wide enough to include all the frequencies we expect to observe due to 
the Doppler shift from the given source, during the time of observation. 

• After the bandwidth has been selected, the data (still in the frequency 
domain) should be windowed, to avoid edge effects in the transformed 
data. 
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• The selected data are then transformed to return to the time domain. At 
this stage we must remove the window used in the data when constructing 
the FFT data base, by simply dividing the new time domain data by w(t). 
This operation recovers the original time data (subsampled) because the 
only regions where the division might not work are the edges of the data 
stream, where the value of the function w(t) may be zero, depending on 
the kind of window used (a problem which, of course, has been overcome 
by the overlapping of the FFTs) . 

• If an FFT under consideration is vetoed or if it is missing, then the data 
are set to zero. 

• Each new group of time domain data is appended to the previous groups, 
after elimination of the overlapped data. Since the overlapping concerns 
half the data we eliminate 1/4 of the data at the beginning and end of 
each stream. The first 1/4 data in the first FFT and the last 1/4 in the 
last FFT can be discarded. The data of missing or vetoed periods, set to 
zero as explained above, are appended to the others in the same way. 

• At this stage we have a subsampled time domain data stream, which 
represents the analytical signal associated with the original data. 

• Now we can take into account the Doppler shift and correct the data as 
previously explained. 

• The final step is to calculate the power spectrum from this subsampled 
time domain data (after data windowing in the time domain) 



APPENDIX D: Uncertainty in the source 
parameters 

Uncertainty in the source position parameters 

In this analysis we used the coordinates a = 17.7 hour, 5 = —29.0 deg to define 
the GC. In order to calculate the region of the sky effectively covered by this 
definition it was necessary to study the effect on the analysis of a source not 
being "exactly" in the GC, since the frequency modulation depends on the 
precise location of the source. To get an idea of the problem in 1991 we plotted 
(Fig. pi) the differ ence in the observed frequencies on Earth between a signal 



from the GC and signals coming from sources at nearby coordinates. 

From the graph it is easy to see that a difference in the right ascension 
of ±0.01 hours leads to a maximum difference in the observed frequency of 
~ 2 • 10~ 4 Hz. This mismatch is at a maximum twice a year, at the beginning 
of June and at the beginning of December. Thus, we studied the effect of the 
mismatch during a run in December, when it was at a maximum. 
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Figure 10: The graph shows the difference in the observed frequency on Earth between 
a signal in the GC and signals coming from nearby coordinates. The x-axis are days 
of 1991. The y-axis are the frequencies, from —2.5 • 10~ 4 to 2.5 • 10~ 4 Hz 



Table I shows the results for uncertainties both in right ascension and in 
declination. The first column gives the error in right ascension (a in hours); 
the second, the error in declination (5 in degrees); the third, the energy of the 
spectrum of the signal, 1/Hz, in the frequency bin of its maximum and in the 
previous and next bins nearest to the maximum; the fourth column gives the 
difference in the frequency of the signal compared to the nominal, expressed 
in number of bins (one bin is 8.1 • 10~ 7 Hz). Fig. 11 is the corresponding 3- 
dimensional plot. The z-axis is the energy of the signal, integrated over the 
three bins. 
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Figure 11: 3-dimensional plot of the data in Table I. The z-axis is the energy of the 
signal integrated over the three bins (maximum,previous, next). The x and y-axes 
are the mismatch in right ascension (hours) and in declination (degrees) 



It is not easy to arrive at a general conclusion, because the final effect 
depends on the mismatch on both the parameters and, in some cases, when the 
two parameters act in opposite direction, the final result is better compared to a 
mismatch in only one of the two parameters. This is why, for example, the result 
when Aa=0.01 hour and A<5=0.3 deg is better than Aa=0.01 hour and A5=0 
deg. Anyway, assuming the analysis is valid even when there is a worsening by 
a factor of 2 in SNR we may conclude that the region of the sky under study is 
definitely within either the volume a = 17.7 ± 0.01 hours and S = —29.0 ± 0.05 
degrees (0.01 hours=0.15 degrees), or the volume a = 17.7 ± 0.005 hours and 
8 = —29.0 ± 0.2 degrees We also note that the observed frequency shift of 
the maximum is consistent with that which would be expected (compare, for 
example, the results in Table I with those in Fig. 10). 
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Table 1: Results for uncertainties in right ascension and in declination 
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Table 2: Results for uncertainty in the source frequency. The first column gives the 
error in the frequency (Au in Hz); ; the second, the energy of the spectrum of the 
signal, 1/Hz, in the frequency bin of its maximum and in the previous and next bins 
nearest to the maximum; the third column gives the difference in the frequency of 
the signal compared to the nominal, expressed in number of bins (one bin is 8.1 ■ 10~ 7 
Hz) 

Table II 



Ay signal energy Af (n bin) 

0.0 (0.46) 1.00 (0.17) +0 

+0.1 (0.29) 0.92 (0.39) +1 

-0.1 (0.71) 0.73 (0.15) -1 



Uncertainty in the source frequency 

To test the extent to which the analysis depends on knowledge of the intrinsic 
frequency of the source, we did a simulation of a spectrum of 14.1 days, by 
introducing a signal at 921.3 Hz and correcting it, during Doppler removal, 
using 921.2 and 921.4 Hz (that is with an error of ±0.1 Hz). From Table II it 
is easy to see that there are no significant differences in the resulting spectra. 
Thus the final result is only slightly affected by even a very "big" error such 
as this. This finding is important because it allows us to analyze just a set of 
discrete frequencies, for example just 1/100 of the frequencies in the original 
FFTs (0.419... mHz). It can be shown that this property is intrinsic to the 

nature of the Doppler correcting factor (eq. |n]), which depends on the difference 
between the intrinsic and the observed frequency. Fig.|l^ shows, for example, 
the difference (over two years) between the terms u\ — and v% — z^2cZ> where 
v\ = 921.300 Hz, V2 = 921.257 Hz, v\& and vid are the observed frequencies 
at the detector, assuming the source to be in the GC. This difference is of the 
order of 3 • 10~ 5 . 
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Figure 12: The top figure shows the observed v xd frequency-921 Hz, assuming a source 
in the GC emitting at z/i=921.300 Hz. The x-axis is in days. The second figure shows 
the same for a signal at z/2=921.257 Hz. The third and the forth are v\ — v\d and 
v 2 — v 2d- The bottom figure is the difference between the third and the fourth figures. 
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